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We investigate the influence of quantum effects on the dielectric and piezoelectric 
properties of barium titanate in its (low-temperature) rhombohedral phase, and show 
the strongly anharmonic character of this system even at low temperature. For 
this purpose, we perform path-integral molecular-dynamics simulations under fixed 
pressure and fixed temperature, using an efficient Langevin thermostat-barostat, 
and an effective hamiltonian derived from first-principles calculations. The quantum 
fluctuations are shown to significantly enhance the static dielectric susceptibility (~ 
by a factor 2) and the piezoelectric constants, reflecting the strong anharmonicity 
of this ferroelectric system even at very low temperature. The slow temperature- 
evolution of the dielectric properties observed below « 100 K is attributed (i) to 
zero-point energy contributions and (ii) to harmonic behavior if quantum effects are 
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turned off. 

I. INTRODUCTION 

It is a common result of solid state physics that, below its Debye temperature, #d, a solid 
exhibits a behavior that significantly deviates from the predictions of classical mechanics, 
because of quantum fluctuations associated to atomic motions: in such conditions, one or 
several vibration modes of the system do not behave classically, i.e. the energy quantum 
Hoj separating their eigenstates is significantly higher than the thermal energy /cgT. In har- 
monic, or mildly anharmonic systems, quantum effects can be theoretically treated in the 
framework of the harmonic or quasiharmonic approximation, using the standard quantiza- 
tion of atomic vibrations through the introduction of the normal coordinates, and many 
solids can be considered as harmonic below a certain temperature. Such harmonic systems 
usually exhibit weak dependence of their dielectric properties with temperature, and their 
dielectric permittivity is not sensitive to quantum effects^. 

However, the influence of quantum effects is more complex in systems in which the micro- 
scopic degrees of freedom move inside a strongly anharmonic energy landscape. For example, 
in a number of crystals exhibiting a form of ferroic order related to atomic displacements 
(ferroelectric, ferroelastic, antiferrodistortive, ferrotorroidic, etc), most of the degrees of free- 
dom typically evolve throughout a multiple-well energy surface. If the system is considered 
as classical (quantum effects neglected), there always exists a temperature below which such 
system exhibits harmonic behavior (it is finally trapped in one single well of the energy 
surface), but in a quantum system, if the quantum fluctuations on the displacements in 
the ground state extend beyond the harmonic region associated to each well, the standard 
treatment of quantum effects by means of the harmonic approximation becomes irrelevant. 
In such cases, where quantum fluctuations of the ground state cause the system to probe the 
anharmonic part of the potential, deviations from the harmonic approximation are expected 
down to zero K. 

Ferroelectric systems are typically characterized by complex multiple-well energy land- 
scapes affecting the polar degrees of freedom, and the interplay between quantum fluctua- 
tions and polar modes at low temperature is strongly system-dependent. In the case of a 
deep multiple-well surface, the effects might be limited to weak quantum derealization. But 
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in the case of a shallow multiple-well surface, dramatic consequences can be observed, up to 
strong tunneling effect that can be able to fight against the order parameter, and even make it 
disappear. One of the most spectacular interaction between quantum fluctuations and polar 
degrees of freedom is indeed the suppression of the ferroelectric polarization due to quantum 
zero point motions in the so-called "quantum paraelectric" materials such as SrTiOgF^ and 
KTaOsj^two systems in which the existence of an underlying ferroelectric instability (thus 
associated to a shallow multiple-well energy surface) is experimentally infered from the very 
large and saturated values of the dielectric permittivity just above zero K, and suggested 
from density-functional theory (DFT) calculations in the case of SrTiO^. Similar effects are 
suggested in another family of compounds, the "high-temperature" quantum paraelectrics, 
such as CaTiOgP^, Lai/2Nai/2Ti03 8 ^ or in a more general way REi/2Nai/2Ti03, 10 in which 
RE features a rare-earth element. However, even in "standard" ferroelectric crystals such as 
BaTiC>3 (BTO), the quantum effects have been shown by path-integral Monte-Carlo tech- 
nique to significantly decrease the phase transition temperatures by w 30-50 KP and to 
strongly modify the shape of the pressure-temperature phase diagram^] U p to room temper- 
ature. Quantum effects can thus strongly influence the structural and dielectric properties 
of ferroelectric (FE) systems, not only at low temperature but also at room temperature 
and beyond. 

Anharmonicities in BTO are impacting the physics even at low temperature, as shown 
by its lowest phase transition (rhombohedral-orthorhombic - expt: 183 KP^, ~ 190 K with 
the hamiltonian used in the present worlsP^). In a simple picture, each phase transition in 
BTO corresponds, upon heating, to the temperature at which the local modes (dipoles) get 
out of the potential well(s) in which they were confined, and come to visit new potential 
energy minima, giving rise to a new value and direction of the macroscopic polarization. 
In other words, phase transitions are the points at which anharmonicities in the potential 
energy surface strongly manifest and impact the physics. This simple microscopic picture 
of BTO has been theorized more than 40 years ago through the well-known "eight-site" 
order-disorder model of Comes et a/P^. In this model, the local dipoles evolve among eight 
off-center satellite sites located along the <lll>-type directions. In the paraelectric phase, 
all the sites are visited with equal probability (resulting in a zero polarization), while in 
the FE phases, only a subset of these sites (the same at all cells due to strong intersite 
correlations) is visited, giving rise to a non-zero polarization. Modern calculations using an 
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effective hamiltonian fitted on first-principles calculations^ show that, in the paraelectric 
phase, the local density of probability is in fact quasi-isotropic (with slight maxima along 
the <111> directions) since the dipoles also spend a significant part of their time between 
these off-center site d 16 * 17 *. At any rate, this simple model confirms the strongly anharmonic 
character of this system, at least above its first phase transition and thus, the complete 
impossibility to describe its phase transitions in the harmonic approximation. Let us mention 
that the signature of anharmonicities is seen also below the first phase transition: in the 
rhombohedral phase, order-disorder mechanims associated with local reversal of the dipoles 
have been observed, both experimentally^ and by calculations,^^ suggesting anharmonic 
behavior also in the ground state phase of BTO. 

On the other hand, the Debye temperature of BTO is commonly placed about 150 K above 
room temperature (#d ~ 480 K?^), showing that the dipole dynamics should be impacted 
by quantum effects up to such high temperature. Therefore, since the structural phase 
transitions occur well below #d, subtle interplay between quantum fluctuations and anhar- 
monicities are expected in BTO. Let us also mention the fact that the so-called "temperature 
rescaling" method^ does not apply to BTO since the effective temperature calculated from 
the phonon density of states obtained at T = K in the rhombohedral phase is about 250 K, 
which falls above the two first phase transitions of BTO, namely in the tetragonal phase. 

In the present work, we investigate the anharmonicities of barium titanate at low tem- 
perature (in its rhombohedral phase), that manifest through the quantum fluctuations of 
its ground state. We describe the influence of quantum effects on the order parameter, on 
the dielectric permittivity and on the piezoelectric constants of BTO, by using path-integral 
molecular-dynamics simulations. 

II. COMPUTATIONAL DETAILS AND THEORETICAL BACKGROUND 

A. Path-integral formalism 

The quantum effects related to atomic motions are accounted for by using the path 
integral (PI) formalism. In this formulation of quantum statistical mechanics, the canonical 
partition function Z is written as a discretized imaginary time path integral. For a quantum 
system containing N (discernable) particles of mass m, Z can be expressed according to: 
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Z = Km ( 2vrmPfc B T )3jVP/2 A f ^/^/({^-^Md^d)...^. (1) 

This integral involves P (Trotter number) replicas of the system labeled by the integer s, 
each replica consisting of a set of N positions = (u ( f\..u^). These replicas characterize 
the discretization of the PI in imaginary time (imaginary time slices). (3 is the statistical 
temperature, (3 = The effective potential V e ff{{U^---U^Y) couples the positions of 
the P slices through: 

p N 

v eff (0^...u^}) = \k{^ - 4 s+l) r + ^({#})]- (2) 

S=l 8=1 

The harmonic term of this effective potential involves a spring constant k = mP ^B T _ 
is the physical potential energy computed inside each imaginary time slice s. Each particle 
(i, s) is thus interacting by harmonic forces with the particles (i, s + 1) and (i, s — 1), forming 
a ring that closes onto itself (u- P+1 ^ = up). 

In the limit of infinite Trotter number, Eq. tends to a functional integral (imaginary- 
time path integral): 



D[/e4f[ r (|)+*(t?W)]* ) ( 3 ) 

the integral being over all paths [r e [0; /3h] — > U(t) G R 3N ] with the cyclic condition 
[U(0) = U(/3h)}, T and $ in the Euclidean action being the kinetic and potential energies, 
that respectively depend on the momenta ^p- and positions U(t). 

The discretized expansion of Eq. [I] is at the root of a formal analogy^ ("classical 
isomorphism"^) between any quantum system and an equivalent classical system made 
of P images of the initial set of N particles, because the multidimensional integral of Eq. [l] 
can be viewed as the canonical partition function of this equivalent classical system. Each 
quantum particle is associated to a ring polymer of P classical particles, these classical par- 
ticles interacting with each other through the "true" physical forces (divided by P) inside 
each slice, and through harmonic forces (acting between each particle of slice (s) and the 
corresponding particles of slice (s — 1) and (s + 1)). Each set of harmonic interactions is 
assumed to close onto itself, forming a closed ring (P+l — > 1). In the limit where the Trotter 
number P — > oo, this equivalent classical system has exactly the same partition function as 
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that of the quantum system under study. The extension to the isothermal-isobaric ensemble 
(NPT) is straightforward^. 

As a consequence, classical simulation techniques such as Monte Carlo (MC) or molec- 
ular dynamics (MD) can be applied to the classical equivalent to estimate numerically the 
thermodynamic properties of the quantum system. The corresponding methods are respec- 
tively called path-integral Monte Carlo and path-integral molecular dynamics (PIMD). The 
estimated properties should be converged with the Trotter number, and if this condition is 
fullfilled, it is possible to compute thermodynamic quantities that exactly include all the 
quantum dispersion effects. However, the path-integral expansion of Eq. [T] assumes distin- 
guishible particles. The physical properties computed by the present PI formalism thus do 
not include exchange between particles, which are assumed to obey Boltzmann statistics^. 

Practical application of PIMD raises technical problems, related to the fact that for 
high Trotter number, the forces acting inside the classical equivalent are mainly harmonic. 
In such harmonic system, obtaining ergodic trajectories by using the standard algorithms 
of MD such as Nose-Hoover thermostat is difficult!^. In order to recover ergodicity, two 
different approaches can be employed: 

(i) a deterministic approach using efficient schemes based on thermostat chainP^^, 

(ii) a stochastic approach based on the use of the Langevin thermostapSHZl 

In the present case, we have found very efficient to use the Langevin dynamics, which is 
extremely powerful to produce an ergodic exploration of phase space by introducing at each 
time step a random force that mimics the "noise" observed in the motion of a brownian 
particle. We have successfully tested the scheme on simple systems (ID and 3D harmonic 
oscillator, double well potential, quartic potential) and found an excellent agreement between 
PIMD and the exact result (obtained by analytic formulae or by a numerical solution of the 
Schrodinger equation). 

Another practical difficulty of PIMD is the existence of modes evolving on very different 
time scales^. To circumvent this difficulty, a specific coordinate transformation ("normal 
mode" or "staging") that diagonalizes the harmonic parts of the PIMD forces could be 
employed. An appropriate choice of fictitious masses gives the same time scale for the 
dynamics of each degree of freedom. Instead of using such a coordinate transformation, we 



performed long enough trajectories (see Sec. II D). 
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B. Hamiltonian 

We use the effective hamiltonian of Zhong et alP^, which is derived from first-principles 
density-functional calculations and has been shown to provide an excellent description of 
the thermodynamics of BTO, especially its complex sequence of phase transitions: rhombo- 
hedral(R) - orthorhombic(O) - tetragonal(T) - cubic(C) - and the (first) order of its phase 
transitions. In particular, although the Curie temperature T c is predicted at w 300 K, i.e. 
about 100 K too low with this hamiltonian, it provides a good value for the R-0 phase 
transition temperature (~ 190 K without inclusion of quantum effects^, and ~ 160 K after 
inclusion of quantum effects, to be compared with the experimental value of 183 Kp^j). The 
degrees of freedom of this hamiltonian are the local modes {ui}, the mechanical displacement 
modes {vi} and the (homogeneous) strain tensor {r]i}. Ui is, roughly speaking, the local polar 
displacement inside cell i, related to the local dipolar moment di through an effective charge 
Z*\ di = Z*ui. For simplicity, the mechanical displacement modes - that allow appearance 
of an inhomogeneous component of the total strain and have been shown to be of very weak 
influence in this material - are not accounted for in the present study. The hamiltonian 
$({-Uj}, {r]i}) is thus a function of the local modes and of the (homogeneous) strain 

tensor components {r/i}. It consists of a (local) "onsite" part, a long-range dipole-dipole in- 
teraction term, a term describing short-range interactions between neighboring local modes 
(up to 3 rd neighbor), a (local) term that couples the local mode to the strain and an elastic 
energy^. In the presence of an external electric field E, a term — £\ Z*Ui.E is added. 

In this work, we use in the discussion the mean local mode <u> as order parameter 
(this is a displacement). An unambiguous relationship with the spontaneous polarization P 
can be made by P = Z*^^, where Q is the unit cell volume and Z* the effective charge 
associated with the local mode (9.956 e). The components of <u> are expressed in ao units, 
where ao is the theoretical lattice constant of BTCP^ computed in the framework of the local- 
density approximation to DFT: ao = 7.46 Bohrs. For each temperature, the macroscopic 
order parameter <u> and the homogeneous strain are obtained by averaging over unit cells, 
(real) time steps, and imaginary time slices, allowing to determine the symmetry of the 
phase (R, O, T or C). 
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C. Langevin barostat 

BTO is simulated under fixed (hydrostatic) pressure and fixed temperature conditions. 
Since Langevin dynamics is very efficient to recover ergodicity within the PI formalism in 
the canonical ensemble, we wish to use this method, not only under fixed temperature but 
also under fixed pressure. The extension of the Langevin method to the isothermal-isobaric 
ensemble has precisely been achieved by Quigley and ProbertP^, giving rise to an algorithm 
in which random and friction forces are applied, not only on the atomic coordinates, but 
also on the supercell vectors. We have thus implemented this "Langevin barostat" within 
the PI formalism. In what follows, second-rank tensors are written in bold. The equations 
of motion on local mode % of slice (s) (with mass m) using the Langevin barostat are: 

d lt = m _ 7 ^ + £M _ pom _ J_ Tr (p^ ) (4) 

dt Jl 1 Nf' Wg ^ ' l ' 

with f} s) = -±V A s)$(u[ s \... } ut ] ) - k{T,P){2u { " ] -uf +l) -u { r 1] ) the PIMD force that 
includes the quantum kinetic energy contribution, which takes the form of an harmonic force 
with spring constant k(T,P) = mPk 2 B T 2 /Ti 2 . The term — corresponds to the friction 

~~*(s) 

force of the Langevin thermostat and L\ is the so-called Langevin force, which is randomly 
drawn at each time step in a gaussian of variance v/ ^^f^ ' St being the time step. The 
momentum pf^ is related to the position uf^ by 

while the matrix of the supercell vectors h and its conjugate momentum pg evolve 
according to 

dh p G h (g) 



dt Wg 



and 



^ = V(t)(X - P ext Id) + ±- T ^Id - 7gPg + L G , (7) 
at Nf ^— ' m 

J i,s 

in which V(t) is the supercell volume (that evolves with time), W g is the "mass" associated 
to the barostat, Nf is the number of degrees of freedom, P ext is the external pressure, Id 
is the identity tensor and X is the internal pressure tensor^ES. In the right member of 
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Eq. [7J one recognizes a friction force on the supercell — 7gPg (jG is a friction coefficient 
for the barostat) and a random force Lg, a 3 x 3 matrix whose components are randomly 
drawn at each time step in a gaussian with variance \j 2lGW ^ ^~ This random force on the 
barostat is symmetrized at each time step to avoid global rotation of the supercell during 
the simulation. Let us precise that the same supercell h is common to all the imaginary 
time slices (there is no replication of the supercell). 

The Langevin algorithm is very sensitive to the quality of the random number generator. 
In this study, we have adopted the routines of R. Chandler and P. Northrupl 33 * 34 l Finally, 
the algorithm is very stable and thermalization is fully achieved within « 20000 MD time 
steps. 

D. Molecular dynamics 

The MD simulations are performed using a 12x12x12 supercell with periodic boundary 
conditions. The time step is 5t=1.0xl0~ 15 s. The mass associated to the local mode (not 
important in the classical case for the computation of ensemble averages, but crucial in 
the quantum case) is 39.0 atomic mass units, determined from the force constant matrix 
eigenvector defining the local mode. This value is the same as that used in Ref. The 
external pressure is fixed at -4.8 GPa, as in Ref. [T31 a negative value that compensates the 
underestimation of the lattice constant within the local-density approximation to DFT. The 
Langevin-PIMD equations of motion are integrated within the Verlet algorithm: at each 
step, the new positions at time t + 5t are computed from the positions at t and t — 5t, 
and also from the velocities at t, as required by the equations of motion of the Langevin 
thermostat-barostat. These velocities are calculated self-consistently (a short internal loop 
is performed at each step of the Verlet algorithm) starting from an estimation of the velocity 
taken from Ref. 1361 

Convergence with the number of imaginary time slices has been carefully studied. Our 
tests show that various properties (polarization, strain and phase transition temperatures) 
of the system from T = 120 to 300 K do not exhibit significant change from P = 8 to 
P = 16. Thus the computation of the dielectric and piezoelectric tensors is achieved at low 
temperature (from T = 30K to T = 137K) by maintaining P x T = 120 x 16 = 1920, 
leading to the use of Trotter numbers as large as P = 64 at the lowest temperature studied 
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(T = 30 K). Note that classical mechanics is recovered by setting the Trotter number P to 
1. 

Equilibrium trajectories of typically 3xl0 5 (low T) up to 5xl0 5 (high T) steps are gen- 
erated at each temperature after an equilibration time that can be long at low temperature, 
due to a very slow dynamics. Thus, at low temperature (T < 60 K), an electric field 
along [111] is applied during the equilibration procedure to help the system reach faster its 
equilibrium state. The friction coefficient of the Langevin thermostat is 7 = 0.5 THz. 

Dielectric and piezoelectric tensors have been computed using a finite differences method, 
by directly applying a finite electric field along [001]. Equilibrated trajectories of PIMD steps 
are generated under static electric field and fixed pressure (-4.8 GPa). Different values of 
the electric field, in the range [—5, +5] x 10 6 V/m, are chosen for each temperature. The 
field is weak enough to induce a linear dielectric response, in both the polarization and the 
strain. 

In order to ensure the sampling accuracy, we compare our PIMD averages with the one 
obtained by PIMD under the staging mode transformation. Since the difficulty to obtain 
such a good sampling could occur with high Trotter number, let us consider the most un- 
favourable case (P = 64 at T = 30K). A long NPT trajectory has thus been followed by 
a Rj 100000 step- iVV"T trajectory performed by fixing the strain tensor components to the 
mean values deduced from the NPT trajectory. This second trajectory was computed by 
using both the primitive coordinates, and the staging coordinates. This allows to demon- 
strate that long trajectories, even using primitive coordinates, provide very well converged 
results, as shown in Tab. [TJ It is important to mention that the primitive averages converge 
to the staging ones after « 20000-30000 times steps, whereas our NPT trajectories contain 
more than 300000 time steps. 

III. SPONTANEOUS POLARIZATION 

As a first step, we investigate the phase sequence of BTO by both classical and quan- 
tum simulations in order (i) to evaluate the importance of the quantum contributions as 
a function of temperature, and (ii) to examine whether the low-temperature spontaneous 
polarization is impacted by quantum fluctuations, providing a first idea of low-temperature 
anharmonicities. 
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TABLE I: Averaged local modes and stress at T = 30 K and for P = 64 obtained over (i) a NPT 
trajectory with primitive coordinates, (ii) a NVT trajectory with primitive coordinates, (hi) a 
NVT trajectory with staging coordinates. The strain hxed for the NVT cases is the one averaged 
along the NPT trajectory. 



E x — E y 



E z = 



Ensemble 
Coordinates 
Local mode (ao) u z 



NPT 
prim. 
0.021612 



NVT NVT 
prim. staging 
0.021592 0.021604 



Stress (GPa) o\ 

04 



4.8 
0.0 



4.800383 4.802153 
0.000117 -0.000038 



Er 



E y = 0; E z = 6.4 x 10 6 V/m 



Ensemble 


NPT 


NVT 


NVT 


Coordinates 


prim. 


prim. 


staging 


Local mode (ao) u x 


0.021110 


0.021127 


0.021125 


u z 


0.022722 


0.022729 


0.022731 


Stress (GPa) a x 


4.8 


4.798115 


4.800640 


0"4 


0.0 


-0.000025 


-0.000031 


0"6 


0.0 


-0.000059 


0.000065 



In the classical case (P = 1), the temperature evolution of polarization is displayed on 
Fig. [Tj We find the expected sequence of phase transitions: R - O - T - C, as experimentally 
observed^ and in excellent agreement with the classical Monte carlo calculations of Zhong 
et. alP^ using the same hamiltonian. In particular, the transition temperatures are very 
close to those obtained by these classical Monte Carlo simulations^: the Curie temperature 
is found at f» 295 K, while the R-0 and O-T transition points are obtained at ~ 190 K and ~ 
230 K, respectively. These preliminary calculations illustrate the accuracy of the Langevin 
barostat to sample the (NPT) ensemble. Note that these transition points are subject 
to uncertainties of ~ 10-15 K related to hysteretic phenomena common to the simulation 
of first-order phase transitions. A more accurate determination of the phase transition 
temperatures would require to achieve a finite-size scaling analysis^! which is beyond the 
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scope of the present study. 

In the quantum case, the transition temperatures (see Fig.[T| are lowered compared to the 
classical case by ~ 10-15%, and agree very well to the ones obtain by path-integral Monte 
Carlo by Zhong and VanderbiltP. This behavior is expected since the quantum fluctuations 
destabilize the ferroelectric order. We also observe an important decrease in the spontaneous 
polarization in the three ferroelectric phases (Fig. [T]). At T = 30 K, we find in the quantum 
case (u x ) = (u y ) = (u z ) = 0.0216 a (0.221 C/m 2 ), whereas in the classical case, we have 
0.0257 a (0.263 C/m 2 ). This difference reflects the fact that the mean polar displacement 
in the ground state is not at the minimum of the potential energy. This is due to the 
quantum fluctuations extending to the asymmetric (and thus anharmonic) region of the 
potential energy surface of the system. Indeed, an harmonic system would have its ground 
state symmetric with respect to the minimum energy point. Despite the multidimensional 
character of this energy surface, it is possible to have an idea of its anharmonicities by 
plotting the energy as a function of polarization P, with all the local modes fixed to the same 
value in all the unit cells of the system: ui = u 2 = ... = = u. This simplified static energy 
landscape is displayed in Fig. [2] along three different directions and help understanding the 
effect of quantum fluctuations: since the energy surface is below (resp. above) the harmonic 
approximation along the [111] direction when the polarization is decreased (resp. increased) 
with respect to its value in the energy minimum, the quantum fluctuations decrease the 
value of the spontaneous polarization, as found in the calculations. 

This significant difference between the classical and quantum systems is a signature of 
the importance of anharmonicities in the rhombohedral phase of BTO. In other words, 
quantum fluctuations lower the spontaneous polarization at low temperature by 15-20%. 
This manifestation of quantum mechanics takes place through the strong anharmonicity of 
the potential energy surface of BTO. 

IV. DIELECTRIC PERMITTIVITY 

We now focus on the rhombohedral phase of BTO below 140 K and compute the di- 
electric and piezoelectric tensors in the presence of quantum fluctuations. This calculation 
is achieved, as explained in the computational part, by applying a finite external electric 
field E along the [001] direction. Figjsl shows the temperature evolution of the transverse 
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and longitudinal components of the static dielectric tensor, computed classically (P=l) and 
quantum-mechanically (P x T = 1920). These components are systematically given, and 
discussed, in the rhombohedral reference system (see Appendix). 

First, we note that, independently from the inclusion of quantum fluctuations, the trans- 
verse component, en, is much larger than the longitudinal component, €33. This behavior is 
expected and in good agreement with the density functional perturbation theory (DFPT) 
calculations of Wu, Vanderbilt and HamanrpS, who obtained en = 265 and £33 = 50 at 
T = K. The ratio between these two components can be explained through the curvature 
of the potential energy surface around the rhombohedral minimum. Fig. [2] shows how the 
energy landscape varies as a function of the polarization around this minimum. The [110] 
and [112] directions give insight into en, the [111] direction into €33. The much sharper 
increase observed along the longitudinal [111] direction (Fig. [2]) shows that €33 is smaller 
(smaller polarization fluctuations) than en along the transverse ones. 

We now examine the influence of quantum fluctuations. In both the classical and 
quantum-mechanical cases, our calculations show that the dielectric constants slowly evolve 
with temperature for temperatures below ~ 100 K (Figj3]). This common feature is ex- 
plained differently according to the case. In the classical approach, the system, as explained 
previously, eventually becomes harmonic at sufficiently low temperature, leading to the slow 
temperature-evolution in the dielectric response. In the quantum case, the system rapidly 
reaches its ground state upon cooling, generating a freezing (saturation) of all the physical 
quantities. The slow evolution of the dielectric response at low temperature is therefore 
attributed to the quantum zero-point effects. For higher temperatures (> 100 K) the di- 
electric constants increase and become very large as approaching from below the R-0 phase 
transition. We systematically find that the inclusion of quantum fluctuations enhances the 
dielectric response, approximately by a factor 2 for both en and e 33 . Such differences be- 
tween the quantum and classical descriptions are the signatures of a strong anharmonicity 
in the potential energy landscape of the rhombohedral phase of barium titanate, since in an 
harmonic system, the static dielectric tensor would not depend on the inclusion of quantum 
effects and would be independent on the temperature^. 

Finally, we have also performed measurements of the dielectric constant, eg 3 , on a [001] 
oriented single crystal in the 10-300 K temperature range by using an impedance analyzer. 
These results are given in Fig. [4] and show two anomalies at 170 ± 5 and 260 ±5 K corre- 
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sponding to the R-0 and O-T ferroelectric phase transitions. The PIMD e 33 (see Appendix) 
are shifted in temperature, to take into account the difference between our calculated and 
our experimental R-0 transition temperatures. The agreement with the experiment is good, 
even though the PIMD values are overestimated at low temperatures and underestimated at 
high temperatures in the R phase. Hence, the effective hamiltonian provides a satisfactory 
behavior. 

V. PIEZOELECTRIC COEFFICIENTS 

Figures [5] and [6] display the temperature evolution of the longitudinal, d 33 , tranverse, 
d 3 i, and shear, d 22 and c?24, piezoelectric constants in the R phase of BTO, computed in the 
classical and quantum-mechanical descriptions. The same types of effects are observed as 
for the dielectric constants. 

We note that, independently from the inclusion of quantum fluctuations, du > d 22 > 
^33 > ^31 and the shear piezoelectric constants are one order of magnitude larger than 
the longitudinal and transverse constants. The values obtained at T = 30 K are in good 
agreement with the DFPT calculations of Wu, Vanderbilt and HamanrP^ (performed at 
T = K). Note that the comparison with T = K results, reported in Tab. [TlJ is relevant 
since the properties slowly evolve with temperature below 100 K. 

The components of the piezoelectric tensor have not been experimentally determined in 
the R phase, to the best of our knowledge. Nevertheless, the longitudinal constant along 
[001], (i^, was derived from the slope of the strain versus electric field curves measured 
at T = 173 K^Sl 10 K below the experimental R-0 transition temperature). d^ 3 (see 
Appendix) calculated at T = 137 K (« 20 K below the R-0 transition temperature) is 
rather close to this experimental value (see Tab. [lT| ). 

It is important to mention that the longitudinal coefficient 0^3, along [001], is higher 
than d 33 , along [111]. This strong anisotropy is due to the contribution of the super large 
shear piezoelectric constants, d 22 and d 2 4. Hence, it is interesting to calculate the orienta- 
tion dependence of d 33 in order to identify the direction along which its value is the more 
enhanced. Figure [7] shows the orientation dependence of d^ 3 in the (110) plane along which 
the maximum value is obtained. It corresponds to a direction close to [001] and this behavior 
is clearly temperature-independent. 
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Hence, in the R phase of BTO the direction of enhanced piezoelectricity is different from 
the direction of the polarization. This feature was previously studied in giant-piezoelectric 
single crystals such as PMN-PT 41 * 42 ! or PZN-PT^. Such piezoelectric properties in [001] 
oriented crystal has been attributed to the very large value of the g?24 shear coefficient. 
This property apparently specific to morphotropic compounds is thus observed in a simple 
perovskite like BTO. Indeed, in the case of BTO the ^24/^33 ratio is around 21 and close to 
the ratio of 22 found for PMN-PT in its R phase. 

TABLE II: Piezoelectic constants (pC/N) of barium titanate obtained by calculations. 

T(K) d 31 d 33 d 22 d24 ^3 [001] 
present work 30 7.5 11 75 236 137 
DFPTE21 6.8 15 70 243 137 
expP 173 275 
present work 137 315 ± 20 



VI. DISCUSSION 

It is useful to compare qualitatively the case of BTO to other ferroelectric systems. Thus 
we now discuss the importance of quantum fluctuations and their interplay with anhar- 
monicities in such systems. We limit the discussion to standard FE systems, that do not 
exhibit any other order parameter than polarization, and to their lowest-temperature phase. 
Let us denote by E zp the zero-point energy associated to atomic motions, and by AV the 
typical (free) energy barrier to overcome to reverse the polarization (as would be given for 
instance by a phenomenological treatment within Landau theory). Different prototypical 
cases can be described. 

(i) FE crystals with very deep double-well free energy landscape exhibit large barrier 
height (AV) compared to E zp and thus a Curie temperature T c relatively high. The typical 
example of such systems is PbTiOa. Upon cooling, for T << T c , the crystal eventually 
behaves as an harmonic system. Therefore, its dielectric permittivity is expected to be 
rather flat at low temperature, and not sensitive to quantum effects. 
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(ii) At the opposite, if E zp > AV or if AV and E zp have similar values, the ground 
state is expected to extend over the different energy minima. In such cases, the ferroelectric 
transition might be suppressed by quantum zero-point fluctuations and the system remains 
paraelectric for all temperatures, despite the existence of a FE instability. This is typically 
the case of the quantum paraelectric crystals, such as KTaOsPln such systems, the dielectric 
permittivity increases upon cooling but, since the phase transition does not take place, it 
eventually saturates to a large value down to zero K. 

(hi) An intermediate situation corresponds to a zero-point energy E zp < AV, but large 
enough anyway so that the ground state extends up to the anharmonic region of the poten- 
tial (while staying confined in a single energy minimum). In such a case, all the physical 
quantities exhibit therefore anharmonic behavior down to T = K. The dielectric permit- 
tivity, in particular, saturates at low temperature, as a property of the ground state. In 
the present work we have demonstrated that BTO belongs to this last case, and is thus an 
anharmonic system down to zero Kelvin. It is important to point out that in the case of the 
classical treatment, harmonic behavior is observed at sufficiently low T, leading here again 
to a plateau in the low-temperature evolution of the dielectric permittivity. However, the 
difference of origin between the quantum and classical plateaus is reflected by their different 
values (the quantum plateau is higher than the classical one). 

VII. CONCLUSION 

In this work, the PIMD simulations have been performed to study low-temperature (R 
phase) dielectric and piezoelectric properties of BTO: spontaneous polarization, dielectric 
susceptibility and piezoelectric constants. We have shown that these three properties are 
different whether they are treated classically or quantum-mechanically. More precisely, 
significant enhancement of the dielectric tensor components and of the piezoelectric constants 
are observed as a consequence of the inclusion of quantum fluctuations. Such enhancement 
is attributed to the anharmonic contributions to the ground state, in which the system 
saturates at low temperature. 

By contrast, without including the quantum effects, the system eventually becomes har- 
monic at low temperature. In BTO, the polarization quantum-fluctuations in the ground 
state therefore extend over a region in which the potential energy surface is strongly an- 
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harmonic. This is corroborated by the significant difference between the low-temperature 
spontaneous polarizations computed classically and quantum-mechanically. 

BTO is a strongly anharmonic system down to zero Kelvin, the anharmonicity being re- 
tained at low temperature by the quantum zero-point effects. The anharmonicity should thus 
be accounted for to achieve realistic predictions in this system, even in its low-temperature 
rhombohedral phase. 



Appendix A: Piezoelectric-coefficient calculations 

We choose the electric field E and the stresses a as independent variables. The electrome- 
chanical properties are therefore described by the dielectric constants at constant stress e?-, 
and the piezoelectric constants d{ a . Since there is no ambiguity, the superscripts a will be 
omitted for simplicity in the following. 

The dielectric and piezoelectric matrix, and di a , of the rhombohedral single domain 
state using the [110], [112] and [111] axis and according to the R3m symmetry is as below 



en 
en 
e 33 

d 24 -2rf 22 

-d 2 2 d 2 2 d 2 4 

d 31 d 3 i d 33 

The coefficients e*j and d* a obtained by axis rotation, i. e. in the referential [100], [010] 
and [001], are as below 
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where 



4i = I fas - en) (Al) 



3 

633 = I (^33 + 2eu) (A2) 



; " ; V3 



(-V2c? 22 - d 24 + 2rf 3 i + ^33) (A3) 
rf* 3 = — (2V2d 2 i + 2rf 24 + 2d 3 i + 4s) (A4) 
d* u = -j= (-2V2d 22 + d 2A - 2d 31 + 2rf 33 ) (A5) 

d* m = (lV2d 22 - 2d 2A - 2d 31 + 2rf 33 ) (A6) 

We consider, at a given temperature, the rhombohedral phase with a spontaneous po- 
larization, P° where P® = P 2 = P® and spontaneous homogenous strains, 77° where 
Vi = V2 = V3 an d V4 = V5 = Ve- When an electric field E is applied along the [001] 
direction, the polarization and homogenous strain variations are given by 

AP t = el E (A7) 
A Va = d* 3a E (A8) 
From equilibrium trajectories for different values of the electric field, it is easy to derive 



the six dielectric and piezoelectric constants, e,,- and di a , by using equations (Al A8) 
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FIG. 1: Temperature evolution of the mean local mode components (proportional to the polariza- 
tion) as obtained by standard MD and PIMD. 
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FIG. 2: Energy landscape of BTO in its rhombohedral phase as a function of the polarization 
variation around the energy minimum, along three directions. 




FIG. 3: Temperature evolution of the longitudinal and transverse components of the dielectric 
tensor, computed using either standard MD or PIMD. The electronic contribution e°° (= 5.24 in 
cubic BTCP^) is not included. 
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FIG. 4: Temperature evolution of the longitudinal dielectric constant 633 along [001]. 
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FIG. 5: Temperature evolution of the c^i and piezoelectric coefficients, computed using stan- 
dard MD and by including the quantum effects through PIMD. 
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FIG. 6: Temperature evolution of the d,22 and c?24 piezoelectric coefficients, computed using stan- 
dard MD and by including the quantum effects through PIMD. 




